function my_simulation(fs,t_initial,n_sample)
h = 1/fs;                   % integration steps
t_final = t_initial + h*n_sample
t=t_initial:h:t_final;		% setup time from 0 to t_final
f=50;						% input frequency
f1=100;                     % input frequency
f2=200;
f3=300;
f4=400;
f5=500;
f6=600;
f7=700;
f8=800;
f9=900;
f10=1000;


n_int = length(t);          % required integration
y(1) = 0;                   % Initial value
dydt(1)=2*pi*f*cos(2*pi*f*t(1)); % Initial acceleration value

% dydt(1)=2*pi*f*cos(2*pi*f*t(1))+ 2*pi*f1*cos(2*pi*f1*t(1)) + 2*pi*f2*cos(2*pi*f2*t(1))...
%     + 2*pi*f3*cos(2*pi*f3*t(1)) + 2*pi*f4*cos(2*pi*f4*t(1)) + 2*pi*f5*cos(2*pi*f5*t(1))...
%     + 2*pi*f6*cos(2*pi*f6*t(1)) + 2*pi*f7*cos(2*pi*f7*t(1)) + 2*pi*f8*cos(2*pi*f8*t(1))...
%     + 2*pi*f9*cos(2*pi*f9*t(1)) + 2*pi*f10*cos(2*pi*f10*t(1));   % Initial acceleration value


for i=2:n_int
    %compute integration
    dydt(2)=2*pi*f*cos(2*pi*f*t(i));
%     dydt(2)=2*pi*f*cos(2*pi*f*t(i)) + 2*pi*f1*cos(2*pi*f1*t(i)) + 2*pi*f2*cos(2*pi*f2*t(i))...
%     + 2*pi*f3*cos(2*pi*f3*t(i)) + 2*pi*f4*cos(2*pi*f4*t(i)) + 2*pi*f5*cos(2*pi*f5*t(i))...
%     + 2*pi*f6*cos(2*pi*f6*t(i)) + 2*pi*f7*cos(2*pi*f7*t(i)) + 2*pi*f8*cos(2*pi*f8*t(i))...
%     + 2*pi*f9*cos(2*pi*f9*t(i)) + 2*pi*f10*cos(2*pi*f10*t(i));   % Input signal
    y(i) = y(i-1) + (h/2)*(dydt(2)+dydt(1));
    %save value to next iteration
    dydt(1)=dydt(2);
end
%show value to console
y(n_int)
figure;
plot(t,y);
axis([0 0.2 -2 2]);
title('Signal with Trapezoidal Integration')
xlabel('time (seconds)')

%compute continuous matlab integration
syms x;
fx = 2*pi*f*cos(2*pi*f*x) ;
% fx = 2*pi*f*cos(2*pi*f*x) + 2*pi*f1*cos(2*pi*f1*x) + 2*pi*f2*cos(2*pi*f2*x)...
%     + 2*pi*f3*cos(2*pi*f3*x) + 2*pi*f4*cos(2*pi*f4*x) + 2*pi*f5*cos(2*pi*f5*x)...
%     + 2*pi*f6*cos(2*pi*f6*x) + 2*pi*f7*cos(2*pi*f7*x) + 2*pi*f8*cos(2*pi*f8*x)...
%     + 2*pi*f9*cos(2*pi*f9*x) + 2*pi*f10*cos(2*pi*f10*x);
expected_value = double(int(fx,x,0,t_final))

err = expected_value - y(n_int)
end